(*  :Title:	Rewrite Rules for Signal Processing Expressions  *)

(*  :Authors:	Brian Evans, James McClellan  *)

(*  :Summary:	To provide system rewrite rules a la Myers and Covell.  *)

(*  :Context:	SignalProcessing`ObjectOriented`RewriteRules`  *)

(*  :PackageVersion:  2.7	*)

(*
    :Copyright:	Copyright 1989-1992 by Brian L. Evans
		Georgia Tech Research Corporation

	Permission to use, copy, modify, and distribute this software
	and its documentation for any purpose and without fee is
	hereby granted, provided that the above copyright notice
	appear in all copies and that both that copyright notice and
	this permission notice appear in supporting documentation,
	and that the name of the Georgia Tech Research Corporation,
	Georgia Tech, or Georgia Institute of Technology not be used
	in advertising or publicity pertaining to distribution of the
	software without specific, written prior permission.  Georgia
	Tech makes no representations about the suitability of this
	software for any purpose.  It is provided "as is" without
	express or implied warranty.
 *)

(*  :History:	*)

(*  :Keywords:	*)

(*
    :Source:	R. H. Bamberger.  {The Directional Filter Bank:
		  A Multirate Filter Bank for the Directional Decomposition
		  of Images},  Georgia Tech Ph. D. Thesis.  November 1990.

		M. M. Covell.  {An Algorithm Design Environment for
		  Signal Processing}.  M.I.T. Ph. D. Thesis.  December,
		  1989.

		R. Crochiere and L. Rabiner, {Multirate Digital Signal
		  Processing}.  Prentice Hall: Englewood Cliffs, NJ.  1983.
 
		B. L. Evans, R. H. Bamberger, and J. H. McClellan,
		  "Rules for Multidimensional Multirate Structures,"
		  to be published in the {IEEE Journal on Signal Processing}.

		J. Kovacevic and M. Veterli, "Perfect Reconstruction Filter
		  Banks with Rational Sampling Rates in One- and
		  Two-dimensions," in Proc. SPIE Conference on Visual
		  Communications and Image Processing (Philadelphia, PA),
		  pp. 1258-1268, November 1989.

		C. S. Myers.  {Signal Representation for Symbolic and
		  Numeric Processing}.  M.I.T. Ph. D. Thesis.  August,
		  1986.  Appendix D.

		P. P. Vaidyanthan and T. Chen, "Some Fundamental Issues
		  in Multidimensional Multirate Signal Processing,"
		  to be published in {IEEE Journal on Signal Processing}.
 *)

(*  :Warning:	*)

(*  :Mathematica Version:  1.2 or 2.0  *)

(*  :Limitation:  *)

(*
    :Discussion:  Systems are defined in the same format as are mathematical
		  operators.  This format was chosen so that Mathematica's
		  powerful pattern matcher could be used without adjustment.
		  A full blown object representation was not needed since
		  no cacheing, fetching, dynamic properties, etc. are needed
		  as they are in signals.  That is, a z-transform operator
		  always exhibits linearity, and its properties do not
		  change.  In addition, systems are not instantiated as
		  are objects.

		  Properties are attached to the head of a system name.
		  The system rewrite rules need these properties to
		  rewrite expressions.

	Imag[x_ y_] :> Real[x] Imag[y] + Imag[x] Real[y]
	Real[x_ y_] :> Real[x] Real[y] - Imag[x] Imag[y]
 *)

(*  :Functions:	 *)



(*  B E G I N     P A C K A G E  *)

BeginPackage[ "SignalProcessing`ObjectOriented`RewriteRules`",
	      "SignalProcessing`ObjectOriented`System`",
	      "SignalProcessing`Support`FilterSupport`",
	      "SignalProcessing`Support`LatticeTheory`",
	      "SignalProcessing`Support`SigProc`",
	      "SignalProcessing`Support`SupCode`" ]


If [ TrueQ[ $VersionNumber >= 2.0 ],
     Off[ General::spell ];
     Off[ General::spell1 ] ];


(*  U S A G E     I N F O R M A T I O N  *)

SPRecursiveRewrite::usage =
	"SPRecursiveRewrite[expr] only returns a fully rewritten version of \
	the signal processing expression expr. \
	The Rewrite knowledge base is let loose on the expression expr \
	without any guidance. \
	The user can see the intermediate expressions by setting the \
	Dialogue option to True or False."

SystemRewriteRules::usage =
	"SystemRewriteRules contains a list of over 50 rewrite rules for \
	operators (systems) . \
	This rule is composed of three other rule bases: \
	(1) rules based on operator properties, \
	(2) rules for single rate digital signal processing, and \
	(3) rules for multirate digital signal processing. \
	The rules in (1) are patterned after those found in Myers' Signal \
	Processing Language and Interactive Computing Environment (SPLICE) \
	and Covell's Algorithm Design Environment (ADE)."

(*  E N D     U S A G E     I N F O R M A T I O N  *)


Begin[ "`Private`" ]


(* commutativemult --  multiplication always commutes in 1-D *)
commutativemult[x_, y_, n_Symbol] := True
commutativemult[x_, y_, n_List] := ( matmult[x, y, n] == matmult[y, x, n] )

(* extractit -- undoes nesting operation *)
extractit[ a_, b__, h_ ] := ( extractit[a, h], extractit[b, h] )
extractit[ h_[a__], h_ ] := ( extractit[a, h] )
extractit[ h_[a__], hother_ ] := h[a] /; ! SameQ[h, hother]

(* gcd  *)
gcd[{m_Integer, n_Integer}] := GCD[m, n]
gcd[m_Integer, n_Integer] := GCD[m, n]

(* getdiagonal *)
getdiagonal[mat_] := getdiagonal[mat, Min[Dimensions[mat]]]
getdiagonal[mat_, n_] :=
	Block [	{i},
		Table[ mat[[i,i]], {i, 1, n} ] ]

(* isDownsamplerQ *)
isDownsamplerQ[ PolyphaseDownsample ] = True
isDownsamplerQ[ Downsample ] = True
isDownsamplerQ[ x_ ] := False

(* isLSIfilterQ *)
isLSIfilterQ[ FIR ] = True
isLSIfilterQ[ IIR ] = True
isLSIfilterQ[ x_ ] := False

(* isUpsamplerQ *)
isUpsamplerQ[ PolyphaseUpsample ] = True
isUpsamplerQ[ Upsample ] = True
isUpsamplerQ[ x_ ] := False

(* iszero  *)
iszero[ l_Symbol ] := SameQ[l, 0]
iszero[ l_List ] := SameQ[l, Apply[Table, {0, {Length[l]}}]]

(* lsi *)
lsi[ x_ ] := LINEAR[x] && SHIFTINVARIANT[x]

(* matmult *)
matmult[a_, b_, n_List] := a . b
matmult[a_, b_, n_Symbol] := a b

(* matinv  *)
matinv[a_, n_List] := Inverse[a]
matinv[a_, n_Symbol] := 1/a

(* mod *)
mod[l_Integer, L_Integer] := Mod[l, L]
mod[l_?IntegerVectorQ, L_?ResamplingMatrixQ] := ResamplingMatrixMod[l, L]

(* modulatorq *)
modulatorq[f_, n_Symbol] := ! FreeQ[f, n]
modulatorq[f_, {n_Symbol}] := ! FreeQ[f, n]
modulatorq[f_, {n_Symbol, nrest__}] :=
	(! FreeQ[f, n]) || modulatorq[f, {nrest}]

(* removecf --  remove common factors in separable or 1-D resampling *)
removecf[m1_, m2_, n_Symbol] :=
	Block [	{gcdvalue},
		gcdvalue = gcd[m1, m2];
		{m1 / gcdvalue, m2 / gcdvalue} ]

removecf[diagmat1_, diagmat2_, n_List] :=
	Block [	{diag1, diag2, gcdlist},
		diag1 = getdiagonal[diagmat1];
		diag2 = getdiagonal[diagmat2];
		gcdlist = Map[gcd, Transpose[{diag1, diag2}]];
		{ DiagonalMatrix[diag1 / gcdlist],
		  DiagonalMatrix[diag2 / gcdlist] } ]

(* resampledfilter *)
resampledfilter[head_, factor_Integer, n_Symbol, coeffs_List, rest___] :=
	Block [	{i, numcoeffs},
		numcoeffs = Length[coeffs];
		newcoeffs = Table[ coeffs[[i]], {i, 1, numcoeffs, factor} ];
		head[n, newcoeffs, rest] ]

(* resampledcoeffs *)
resampledcoeffs[factor_Integer, coeffs_List] :=
	Block [	{i, matchlist, numcoeffs},
		numcoeffs = Length[coeffs];
		matchlist = Table[0, {numcoeffs}];
		For [ i = 1, i <= numcoeffs, i += factor, matchlist[[i]] = _ ];
		MatchQ[coeffs, matchlist] ]

(* sepresq -- detect separable or 1-D resampling *)
sepresq[mat_, n_Symbol] := IntegerQ[mat]
sepresq[mat_, n_List] := DiagonalMatrixQ[mat] && ResamplingMatrix[mat]

(* shiftq *)
shiftq[expr_] := MatchQ[expr, Shift[l_, n_][x_]]

(* switchtest *)
switchtest[L_, M_, Right] :=
	Block [	{cond, dL, dM, uLinv, uMinv, vLinv, vMinv},
		{cond, {uMinv, dM, vMinv}, {uLinv, dL, vLinv}} =
		  SmithFormSameV[SmithNormalForm, M, L];
		cond && RelativelyPrime[dL, dM] ] /;
	ResamplingMatrix[L] && ResamplingMatrix[M]

switchtest[L_, M_, Left] :=
	First[ SmithFormSameU[SmithNormalForm, M, L] ] /;
	ResamplingMatrix[L] && ResamplingMatrix[M]


(*  G L O B A L S  *)

RulesBasedOnProperties = {

	(*  Shift invariant operator applied to shift  [Myers, 216]	*)
	t_ [ Shift[l_, n_][x_] ] :>
		Shift[l,n] [ t [ x ] ] /; SHIFTINVARIANT[t] ,
	t_ [ s1_, s__ ] :>
		Block [	{signals, signalstomove},
			signals = {s1, s};
			signalstomove = Select[signals, shiftq];
			MoveOperatorsToFront[t, signals,
					     signalstomove, Shift] ] /;
		SHIFTINVARIANT[t] && Count[{s1, s}, Shift[l_,n_][x_]] > 0 ,

	(*  Memoryless operator applied to shifts  [Myers, 216]		*)
	(*    t_ [ Shift[l_, n_] [x_] ] :> Shift[l, n] [ t [ x ] ] /;	*)
	(*    MEMORYLESS[t],						*)
	(*  is handled by the one-input to a shift-invariant system	*)
	(*  rule above because one-input memoryless systems are SI.	*) 
	t_ [ Shift[l_, n_] [x_], y__ ] :> Shift[l, n] [ t [ x, GetArgs[y] ] ] /;
		MEMORYLESS[t] && SameFormQ[Shift[l,n][s_], y] ,

	(*  Memoryless operator applied to downsample  [Myers, 216]	*)
	t_ [ Downsample[m_, n_] [ x_ ] ] :> Downsample[m,n] [ t[x] ] /;
		MEMORYLESS[t] ,
	t_ [ Downsample[m_, n_] [ x_ ], y__ ] :>
		Downsample[m,n] [ t [ x, GetArgs[y] ] ] /;
		MEMORYLESS[t] && SameFormQ[Downsample[m, n][s__], y] ,

	(*  Memoryless operator applied to upsample  [Myers, 216]	*)
	t_ [ Upsample[l_, n_] [ x_ ] ] :>
		Upsample[l, n] [ t[ x ] ] /;
		MEMORYLESS[t] ,
	t_ [ Upsample[l_, n_] [ x_ ], y__ ] :>
		Upsample[l, n] [ t [ x, GetArgs[y] ] ] /;
		MEMORYLESS[t] && SameFormQ[Upsample[l, n][s__], y] ,

	(*  Additive system and sum input  [Myers, 217]		*)
	t_ [ x_ + y_ ] :>  t[x] + t[y] /; ADDITIVE[t] ,

	(*  Additive system and zero input  [Myers, 217]	*)
	t_ [ zero_ ] :> 0 /; ADDITIVE[t] && ZeroQ[zero] ,

	(*  Homogeneous system and scaled input  [Myers, 217]	*)
	t_ [ g_ x_. ] :> g t[x] /;
		HOMOGENEOUS[t] && (! SameQ[g, x, 1]) &&
		Implies[ ! AtomQ[t],
			 MyFreeQ[g, GetOperatorVariables[t]] ] ,

	(*  Commutative system reorder  [Myers, 216]	*)
	t_ [ x_, y_] :> t [ y, x ] /; COMMUTATIVE[t] ,

	(*  Associative system reorder  [Myers, 217]	*)
	t_ [ x_, t_ [ y_, w_ ] ] :> t [ t [ x, y ], w ] /; ASSOCIATIVE[t] ,

	(*  Linear shift-invariant props not in CM's thesis [O & W, 85-86]  *)
	t_ [ u_ [ x_ ] ] :> u [ t[x] ] /; lsi[t] && lsi[u]
	  
}


SingleRateRules = {

	(*  Collect cascaded convolution operations under one operator	*)
	Convolve[n_] [ x___, Convolve[n_][y1_, yrest__], z___  ] :>
		Convolve[n] [ extractit[ x, y1, yrest, z, Convolve[n] ] ] ,

	(*  Convolution of upsampled signals [Myers, 221]  *)
	Convolve[n_] [ Upsample[l_, n_][x_], y__ ] :>
		Upsample[l, n] [ Convolve[n] [ x, GetArgs[y] ] ] /;
		SameFormQ[ Upsample[l,n][s__], y ] ,

	(*  Combining redundant operations [Myers, 217/220] [Bamberge]  *)
	Periodic[m1_, w_] [ Periodic[m2_, w_] [x_] ] :>
		Periodic[m1 m2, w][x] ,
	ScaleAxis[l1_, w_] [ ScaleAxis[l2_, w_] [x_] ] :>
		ScaleAxis[l1 l2, w][x] ,
	Shift[m1_, n_] [ Shift[m2_, n_] [x_] ] :>
		Shift[m1 + m2, n][x] ,

	(*  1-D Shift of a periodic signal-- could repeat ad infinitum  *)
	Shift[m_, n_Symbol] [ Periodic[k_, n_Symbol] [x_] ] :>
		Shift[ Mod[m,k], n] [ Periodic [k,n] [x] ] /;
		( ! SameQ[ m, Mod[m, k] ] ) && ( ! MatchQ[ m, Mod[f_, k] ] ) ,

	(*  Conjugate of a product [Myers, 218]  *)
	Conjugate[ x_ y_ ] :> Conjugate[x] Conjugate[y] ,

	(*  Real part of a product [Myers, 218] *)
	Re[ x_ y_ ] :> Re[x] Re[y] - Im[x] Im[y] ,

	(*  Imaginary part of a product [Myers, 218] *)
	Im[ x_ y_ ] :> Re[x] Im[y] + Im[x] Re[y] ,

	(*  Extended math functions [Myers, 220]  *)
	(*  signal + 0 = signal				built-in	*)
	(*  a^n1 a^n2 = a^(n1 + n2)			built-in	*)

	(*  Factorizations [Myers, 221] *)
	x_ + g_ x_ :> (1 + g) x ,
	g1_ x_ + g2_ x_ :> (g1 + g2) x	(* also covers g x + g y = g (x + y) *)
}


MultirateRules = {

	(*  Identities for downsamplers [Myers, 217] [Evans, fig. 6] *)

	Shift[n0_, n_] [ Downsample[M_, n_] [x_] ] :>
		Downsample[M, n] [ Shift[matmult[M, n0, n], n][x] ] ,

	Downsample[M_, n_] [ Shift[n0_, n_] [x_] ] :>
		Shift[matmult[matinv[M, n], n0, n], n] [ Downsample[M, n][x] ] ,

	phi_ Downsample[M_, n_][x_] :>
		Downsample[M, n][x Upsample[M, n][phi]] /;
		modulatorq[phi, n] ,

	Downsample[M_, n_][ phi_ x_ ] :>
		Downsample[M, n][x] Downsample[M, n][phi] /;
		modulatorq[phi, n] ,


	(*  Identities for upsamplers [Myers, 219] [Evans, fig. 7] *)

	Upsample[L_, n_] [ Shift[n0_, n_][x_] ] :>
		Shift[matmult[L, n0, n], n][ Upsample[L, n][x] ] ,

	Shift[n0_, n_][ Upsample[L_, n_] [ x_] ] :>
		Upsample[L, n][ Shift[matmult[matinv[L, n], n0, n], n][x] ] ,

	phi_ Upsample[L_, n_][x_] :>
		Upsample[L, n][x Downsample[L, n][phi]] /;
		modulatorq[phi, n] ,

	Upsample[L_, n_][ phi_ x_ ] :>
		Upsample[L, n][phi] Upsample[L, n][x] /;
		modulatorq[phi, n] ,


	(*  Identities for cascades of up and downsamplers [Evans, fig. 9] *)

	Downsample[S_, n_] [ Upsample[S_, n_][x_] ] :> x ,

	Upsample[S_, n_] [ Downsample[S_, n_][x_] ] :>
		x ImpulseTrain[S, n] ,

	Downsample[M1_, n_] [ Downsample[M2_, n_] [x_] ] :>
		Downsample[matmult[M1, M2, n], n][x] ,

	Upsample[L1_, n_] [ Upsample[L2_, n_] [x_] ] :>
		Upsample[matmult[L1, L2, n], n][x] ,

	Downsample[M_Integer, n_Symbol] [ Upsample[L_Integer, n_Symbol][x_] ] :>
		Upsample[ L / M, n ][x] /;
		IntegerQ[ L / M ],

	Downsample[M_List, n_List] [ Upsample[L_List, n_List][x_] ] :>
		Upsample[Inverse[M] . L, n][x] /;
		IntegerVectorQ[ Inverse[M] . L ] ,


	(*  Commutativity of cascades of up and downsamplers [Evans, fig. 10] *)
	(*  Phase out last two rules when RelativelyPrime is complete  *)

	Downsample[M1_, n_] [ Downsample[M2_, n_][x_] ] :>
		Downsample[M2, n] [ Downsample[M1, n][x] ] /;
		commutativemult[M1, M2, n] ,  (* always commutes in 1-D *)

	Upsample[L1_, n_] [ Upsample[L2_, n_][x_] ] :>
		Upsample[L2, n] [ Upsample[L1, n][x] ] /;
		commutativemult[L1, L2, n] ,  (* always commutes in 1-D *)

	Downsample[M_, n_] [ Upsample[L_, n_][x_] ] :>
		Upsample[L, n] [ Downsample[M, n] [x] ] /;
		commutativemult[M, L, n] && RelativelyPrime[L, M] ,

	Upsample[L_, n_] [ Downsample[M_, n_][x_] ] :>
		Downsample[M, n] [ Upsample[L, n] [x] ] /;
		commutativemult[L, M, n] && RelativelyPrime[L, M] ,

	Downsample[M_, n_List] [ Upsample[L_, n_][x_] ] :>
		Upsample[L, n] [ Downsample[M, n] [x] ] /;
		CommutableResamplingMatricesQ[M, L] ,

	Upsample[L_, n_List] [ Downsample[M_, n_][x_] ] :>
		Downsample[M, n] [ Upsample[L, n] [x] ] /;
		CommutableResamplingMatricesQ[L, M] ,


	(*  Fundamental Identities Based on Smith Form Dec. [Evans, fig. 11] *)

	Upsample[U_, n_List][x_] :> Downsample[Inverse[U], n][x] /;
		RegUnimodularResMatrixQ[U] ,

	Downsample[U_, n_List][x_] :> Upsample[Inverse[U], n][x] /;
		RegUnimodularResMatrixQ[U] ,

	Downsample[V_, n_List] [ Upsample[U_, n_][x_] ] :>
		Upsample[Inverse[V] . U, n] [x] /;
		RegUnimodularResMatrixQ[U] && RegUnimodularResMatrixQ[V] ,

	Upsample[U_, n_List] [ Downsample[V_, n_][x_] ] :>
		Downsample[V . Inverse[U], n] [x] /;
		RegUnimodularResMatrixQ[U] && RegUnimodularResMatrixQ[V] ,


	(*  Removing redundancy in up/downsamplers cascades [Evans, fig. 12]  *)

	Downsample[M_, n_] [ Upsample[L_, n_][x_] ] :>
		Block [	{newL, newM},
			{newL, newM} = removecf[L, M, n];
			Downsample[newM, n][ Upsample[newL, n][x] ] ] /;
		sepresq[M, n] && sepresq[L, n] && ! RelativelyPrime[M, L] ,

	Upsample[L_, n_List] [ Downsample[M_, n_List][x_] ] :>
		Block [	{cond, dL, dM, uLinv, uMinv, vLinv, vMinv},
			{cond, {uMinv, dM, vMinv}, {uLinv, dL, vLinv}} =
			  SmithFormSameV[SmithNormalForm, M, L];
			newM = Inverse[uMinv] . dM;
			newL = Inverse[uLinv] . dL;
			Upsample[newL, n][ Downsample[newM, n][x] ] ] /;
		First[ SmithFormSameV[SmithNormalForm, M, L] ] ,

	Downsample[M_, n_List] [ Upsample[L_, n_List][x_] ] :>
		Block [	{cond, dL, dM, uLinv, uMinv, vLinv, vMinv},
			{cond, {uLinv, dL, vLinv}, {uMinv, dM, vMinv}} =
			  SmithFormSameV[SmithNormalForm, L, M];
			newL = Inverse[uLinv] . dL;
			newM = Inverse[uMinv] . dM;
			Downsample[newM, n][ Upsample[newL, n][x] ] ] /;
		First[ SmithFormSameU[SmithNormalForm, L, M] ] ,


	(*  Switch operations in a non-commutable cascade [Evans, fig. 13] *)

	Upsample[L_, n_List][ Downsample[M_, n_List][x_] ] :>
		Block [	{cond, dL, dM, uLinv, uMinv, vLinv, vMinv},
			{cond, {uMinv, dM, vMinv}, {uLinv, dL, vLinv}} =
			  SmithFormSameV[SmithNormalForm, M, L];
			newM = dM . uLinv;
			newL = dL . uMinv;
			Downsample[newM, n][ Upsample[newL, n][x] ] ] /;
		switchtest[L, M, Right] ,

	Downsample[M_, n_List][ Upsample[L_, n_List][x_] ] :>
		Block [	{cond, dim, dL, dM, i, uLinv, uMinv, vLinv, vMinv},
			{cond, {uMinv, dM, vMinv}, {uLinv, dL, vLinv}} =
			  SmithFormSameU[SmithNormalForm, M, L];
			dim = Length[n];
			For [ i = 1, i <= dim, i++,
			      dcommon = GCD[ dM[[i,i]], dL[[i,i]] ];
			      dM[[i,i]] /= dcommon;
			      dL[[i,i]] /= dcommon ];
			newM = vLinv . dM;
			newL = vMinv . dL;
			Upsample[newL, n][ Downsample[newM, n][x] ] ] /;
		switchtest[M, L, Left] ,

	(*  Downsampled of shifted upsample [Myers, 219]  *)
	Downsample[L_, n_] [ Shift[l_, n_] [ Upsample[L_, n_][x_] ] ] :>
		Block [ {n0},
			n0 = matmul[matinv[L, n], l, n];
			Shift[n0, n][x] ] /;
		iszero[ mod[l, L] ] ,

	Downsample[L_, n_] [ Shift[l_, n_] [ Upsample[L_, n_][x_] ] ] :>
		0 /;
		! iszero[ mod[l, L] ] ,

	(* Downsample of shifted upsample with co-prime M and L [Myers, 219] *)
	Downsample[M_, n_] [ Shift[l_, n_] [ Upsample[L_, n_][x_] ] ] :>
		Block [	{l1, l2},
			{l1, l2} = EuclidFactors[l, M, L];
			Shift[l1, n] [
			  Upsample[L, n] [
			    Downsample[M, n] [
			      Shift[l2, n] [ x ] ] ] ] ] /;
		RelativelyPrime[M, L, Right] ,

	(*  Interleave as sum of upsampled signals  [Myers, 219]  *)
	Summation[i_, 1, L_, 1] [
	  Shift[i_, n_] [
	    Upsample[L_, n_] [
	      Element[x_List, i_] ] ] ] :>
		Interleave[n][ GetArgs[ToList[x]] ] /;
		SameQ[ Length[x], L ] ,


	(*  1-D Downsample of interleave (L rel. prime to m)  [Myers, 220]  *)
	Downsample[m_, n_Symbol] [Interleave[n_Symbol][x1_, x__]] :>
		Block [	{ii, len, xlist},
			xlist = List[x1, x];
			len = Length[xlist];
			oneterm = Downsample[m, n] [
				    Shift[- ii m / len, n] [
				      Element[ xlist, Mod[ii m, len] + 1 ] ] ];
			Apply[Interleave[n], Table[oneterm, {ii,0,len-1}]] ] /;
		gcd[Length[{x1, x}], m] == 1 ,

	(*  Switching filter and decimator/interpolator [C & R]  *)
	(*  -- pattern matching assumes that the Downsample and  *)
	(*     and PolyphaseDownsample takes the same first and  *)
	(*     last parameter, but the middle parameter(s) vary  *)
	filter_[n_, g_List, rest___][ upsample_[l_, args___, n_][x_] ] :>
		upsample[l, args, n][
			resampledfilter[filter, l, n, g, rest][x] ] /;
		isLSIfilterQ[filter] && isUpsamplerQ[upsample] &&
		resampledcoeffs[l, g] ,

	downsample_[m_, args___, n_][ filter_[n_, g_List, rest___][x_] ] :>
		resampledfilter[filter, m, n, g, rest][
			downsample[m, args, n][x] ] /;
		isLSIfilterQ[filter] && isDownsamplerQ[downsample] &&
		resampledcoeffs[m, g] ,


	(*  Polyphase forms of FIR interpolation and decimation structures  *)
	FIR[n_, coeffs_, rest___][ Upsample[l_, n_][x_] ] :>
		PolyphaseUpsample[l, FIR[n, coeffs, rest], n][x] ,

	Downsample[m_, n_][ FIR[n_, coeffs_, rest___][x_] ] :>
		PolyphaseDownsample[m, FIR[n, coeffs, rest], n][x]

}


(*  SystemRewriteRules --  62 rules (systems only; no signals allowed)  *)
SystemRewriteRules =
	Join[RulesBasedOnProperties, SingleRateRules, MultirateRules]




(*  A P P L Y I N G     T H E     R E W R I T E     R U L E S  *)

(*  The Rewrite rules base can be driven by several routines.  *)

(*  Recursive Rewrite rules or when all else fails	*)
(*	op_[p__][args__] :> op[p] [ f[args] ] ,		*)
(*	op_[p__] :> op[p] ,				*)

Options[ SPRecursiveRewrite ] := { Dialogue -> False }

(*  SPRecursiveRewrite   *)
SPRecursiveRewrite[e_, options___] :=
	rewrite[e, Join[ToList[options], Options[SPRecursiveRewrite]]]

rewrite[e_, oplist_] := rewrite[e, oplist, Replace[Dialogue, oplist]]

rewrite[e_, oplist_, False] := ReplaceRepeated[e, SystemRewriteRules]
rewrite[e_, oplist_, True] := RewriteTracer[e, oplist]
rewrite[e_, oplist_, All] := RewriteTracer[e, oplist]
rewrite[e_, oplist, dialogue_] := Message[Dialogue::notvalid, dialogue]


(*  RewriteTracer  *)
RewriteTracer[e_, oplist_] :=
	Block [	{difference = True, oldexpr, newexpr},
		newexpr = e;
		Print[e];
		While [ difference,
			oldexpr = newexpr;
			newexpr = Replace[oldexpr, SystemRewriteRules];
			difference = ! SameQ[oldexpr, newexpr];
			If [ difference,
			     Print["which becomes"]; Print[newexpr] ] ];
		newexpr ]


(*  E N D     P A C K A G E  *)

End[]
EndPackage[]

If [ TrueQ[ $VersionNumber >= 2.0 ],
     On[ General::spell ];
     On[ General::spell1 ] ];


(*  W R I T E     P R O T E C T I O N  *)

Block [	{newfuns},
	newfuns = { SPRecursiveRewrite };
	Combine[ SPfunctions, newfuns ];
	Apply[ Protect, newfuns ];
	Protect[ SystemRewriteRules ] ]


(*  E N D I N G     M E S S A G E  *)

Print["System rewrite rules have been loaded."]
Null
